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ABSTRACT 


An improved hybrid particle-finite element method has been developed for hypervelocity 
impact simulation. The method combines the general contact-impact capabilities of particle 
codes with the true Lagrangian kinematics of large strain finite element formulations. Un- 
like some alternative schemes which couple Lagrangian finite element models with smooth 
particle hydrodynamics, the present formulation makes no use of slidelines or penalty forces. 
The method has been implemented in a parallel, three dimensional computer code. 

Simulations of three dimensional orbital debris impact problems, using this parallel hybrid 
particle-finite element code, show good agreement with experiment and good speedup in 
parallel computation. The simulations included single and multi-plate shields as well as 
aluminum and composite shielding materials, at an impact velocity of eleven kilometers per 
second. 
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1 Introduction 


The orbital debris hazard to the International Space Station (ISS) and other space structures 
has focused a significant research effort on the problem of spacecraft shielding design. To 
date shield design work has relied primarily on experimental hypervelocity impact research 
(Christiansen et al., 1999), with impact simulation playing a relatively minor role. However 
a number of factors suggest that future design work will place increased emphasis on the use 
of simulation: 

• Experimental studies using light gas guns (LGG) and inhibited shaped charge launchers 
(ISCL) are at present unable to investigate the entire projectile velocity and kinetic 
energy range of interest 

• The increased use of composite materials (in both shielding and aerospace structures) 
and the introduction of multi-layered shields has greatly expanded the number of ex- 
periments required to fully investigate each shielding design problem 

• An emphasis on faster, lower cost design methods in a maturing aerospace industry 
motivates the expanded use of computer aided design tools 

• The increased availability of high performance and massively parallel computing hard- 
ware will expand the range of design problems which can be investigated through 
simulation 

Despite these strong motivations, progress in the application of impact codes to spacecraft 
shielding design problems has been relatively slow. The most important reason is that the 
numerical methods embodied in traditional continuum Lagrangian and Eulerian codes are 
not well suited to address certain noncontinuum physics associated with the shielding design 
problem (Fahrenthold, 1998). Hence initial work on the application of conventional contin- 
uum codes to debris shield impact simulations has been followed by a period of research, 
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during which alternative numerical methods have been tested, modified, and / or developed 
to better address the problem of interest. 

Most recent work simulating orbital debris impact effects has employed either pure parti- 
cle or mixed particle-continuum methods (Hayhurst et al., 1998 and Fahrenthold and Horban, 
1999) since only particle-based kinematic schemes offer both an efficient solution to the de- 
bris propagation problem and an entirely general representation of contact-impact. Work 
based on pure particle methods has encountered difficulties with accurate modeling of mate- 
rial strength effects (Faraud et al., 1999), and other complications (Libersky et al., 1997). It 
appears that some mixed or hybrid particle-continuum method will prove most effective in 
meeting the need for fundamental improvements in simulation-based design of orbital debris 
shielding. 

The present report describes work performed to evaluate a particular new hybrid particle- 
continuum method (Fahrenthold and Horban, 2000), developed to simulate orbital debris 
impact problems. The numerical method is evaluated here, via simulation of a set of ISCL 
experiments, the latter conducted by Grosch (1996 and 1997) to investigate the performance 
of ISS shielding in oblique impacts at a velocity of eleven kilometers per second. The simula- 
tions discussed include: (1) both Whipple and multi-plate shield designs, (2) both aluminum 
and composite shielding materials, and (3) both hollow cylindrical projectiles (produced by 
the ISC launchers) and mass equivalent spherical projectiles (for comparison to lower velocity 
LGG tests). 

The simulations were performed using a parallel code written by the author. In addition 
to the simulation results, speedup data is presented for test problems run on up to 128 pro- 
cessors, on an Origin 2000 system operated by the Numerical Aerospace Simulation facility 
at NASA Ames Research Center. 
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2 Numerical method 


A detailed description of the hybrid particle-finite element technique used here is provided 
in later sections. A brief summary of the numerical method follows. 

In the hybrid formulation employed here, particles and finite elements are used simulta- 
neously but not redundantly to represent different physical effects. The particles are used 
to represent all inertia effects as well as the thermomechanical response of the medium in 
compressed states. The particle center of mass coordinates in the reference configuration de- 
fine Lagrangian finite elements, which are simultaneously employed to represent interparticle 
forces associated with tension and elastic-plastic shear. Damage variables are introduced as 
internal states for the finite elements, and evolve with the material history to represent the 
loss of tensile and shear strength and stiffness under thermomechanical loading. Element 
failure due to spall, melting, accumulated plastic strain or other physical criteria results in 
the loss of interparticle forces associated with element shear and tension, so that particles 
unassociated with any intact elements are free to flow under contact-impact loads. No mass 
or energy is discarded at element failure, and no rezoning is required to model the transition 
from an intact to a fragmented medium. 

This hybrid modeling technique avoids the tensile instabilities and numerical fracture 
problems which arise with some pure particle methods, the use of slidelines and penalty 
forces which characterize pure Lagrangian finite element methods, and the mixed material 
thermodynamics and numerical diffusion which are features of Eulerian techniques. As 
indicated in the sections which follow, this hybrid particle-element methodology can be 
effectively applied in the simulation of rather complex three dimensional debris shielding 
problems. 
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3 Hypervelocity impact simulation 

Most computer simulation work on hypervelocity impact problems has employed either La- 
grangian finite element techniques (Hallquist, 1983) or Eulerian finite difference or finite 
volume techniques (McGlaun et ah, 1990). The advantages and disadvantages of such meth- 
ods are well known, and have led to the application of different codes to distinct problem 
classes. Recent numerical methods research (Belytschko et ah, 1996 and Liu et ah. 2000) has 
included a focus on the development of new Arbitrary Lagrangian-Eulerian (ALE. Budge 
and Peery, 1993), smooth particle hydrodynamic (SPH, Stellingwerf and Wingate, 1993), 
and mixed particle-continuum methods (Hayhust et ah, 1998 and Fahrenthold. 1998) for 
problems which are difficult to address with traditional Lagrangian or Eulerian techniques. 
An example of the latter class of problems is the simulation of hypervelocity impact on 
orbital debris shielding (Christiansen et ah, 1999). 

Although significant advances have been made in the formulation and improvement of 
ALE and SPH techniques, it appears from work to date (Libersky et ah. 1997 and Faraud et 
ah, 1999) that some mixed particle-continuum formulation is best suited to address the or- 
bital debris shielding problem. The most popular type of mixed-particle continuum approach 
is a coupled SPH-finite element formulation. The latter approach has been implemented in 
the EPIC (Johnson et ah, 1993), and AUTODYN (Hayhust et ah, 1998) codes, and in general 
relies on a penalty-based contact-impact algorithm to link conventional SPH and Lagrangian 
finite element models. An alternative mixed particle-continuum formulation was developed 
by Fahrenthold and Horban (1999). The latter approach may be characterized as a hybrid 
(as opposed to a coupled ) formulation, since it makes simultaneous use of both elements and 
particles to model distinct physical effects in all of the impacting materials. 

The hybrid particle-element formulation of Fahrenthold and Horban (1999) employed 
deforming Lagrangian particles and a penalty method to represent contact-impact. In the 
interest of eliminating the penalty treatment of contact-impact effects, Fahrenthold and Koo 
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(2000) modified the formulation, introducing a new kernel based density interpolation for 
compressed states, while retaining an element-based description of tension and shear. They 
demonstrated the effectiveness of this improved hybrid scheme in one dimensional thermoe- 
lastic shock compression and tensile wave propagation problems. The latter formulation was 
extended to the thermoelastic-plastic-damage case by Fahrenthold and Horban (2000). The 
present report outlines the last cited improved hybrid particle-element formulation, intro- 
duces new damage models and a rate dependent yield stress, and evaluates the formulation 
by comparing simulations to published results of hypervelocity impact experiments. 

The principal advantage of the numerical method described here is its seamless integra- 
tion of the general contact-impact capabilities of particle methods with the true Lagrangian 
strength models of finite element formulations. Its principal disadvantage is the computa- 
tional cost of incorporating both element and particle kinematics in a single code, although it 
should be emphasized that nowhere are the particle based and element based computations 
redundant. 

4 Modeling methodology 

Development of the model described here is based on interpolation concepts derived from 
the particle and finite element literature and a Hamiltonian model formulation approach. 
A Hamiltonian methodology replaces the weighted residual solution techniques commonly 
applied in finite element analysis and a variety of (sometimes ad hoc) model formulation 
procedures used in pure particle based modeling. Total entropy variables are introduced here 
as Hamiltonian displacements, so that general thermomechanical dynamics are included in 
the formulation. 

The fixed mass particles used in the present work are non-deforming, and are associated 
with kernel functions used to calculate the continuum density in compressed states. Associ- 
ated with the mass of the ith particle ( m ^) is a total entropy (S^), so that mass density 
(pfi^) and entropy per unit mass (s^) are the thermodynamic states used to calculate the 
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particle pressures and temperatures for any chosen equations of state. 

Here the particle center of mass coordinates (c^) in the reference configuration define 

(j) 

connectivity’s for large strain Lagrangian finite elements, with reference volumes V 0 . In 
the reference configuration the particle packing scheme is body centered cubic, with eight 
corner particles defining a hexahedral finite element, the latter used to describe elastic (E e ^') 
and plastic (E^) strain tensors and thereby quantify elastic-plastic shearing deformation in 
the continuum. Associating the body centered particle in each element with each of the six 
faces of the aforementioned hexahedron defines six five-noded subelements whose reference 
volumes (Vj J,k ^) and Jacobians (J^ ,fc )) are used to quantify hydrostatic tension. 

Normal (D^) and deviatoric (S J ^ ) continuum damage variables are introduced here for 
the finite elements, and are used to quantify reductions in element strength and stiffness, 
as a function of each element’s thermal and mechanical loading history. The damage vari- 
ables evolve from an initial (undamaged) value of zero, to a maximum value of one, the 
latter representing a complete loss of element cohesion. All strain energy release associated 
with damage evolution is accounted for, as irreversible entropy production, and no rezoning 
or mass discard is associated with element failure. Particles unassociated with any intact 
elements are free to flow in response to contact-impact loads. 

The preceding methodology offers a unique combination of features, and is well suited 
to address particular hypervelocity impact applications. The sections which follow discuss 
the particle and element based interpolations, the stored energy functions for the modeled 
system, the dissipative process models, the entropy evolution relations, the final form of 
Hamilton’s equations, and numerical application of the method. 

5 Interpolations 

The density interpolation for compressed states is 

„«)=pM + ^) + pW (1) 
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( 2 ) 


where p^ is the reference density for the ith particle and 
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with ITT.?) and lyTt) kernel functions for reference configuration nearest neighbors and all 
other neighbor particles respectively. For a system of n particles the summation limits satisfy 

n= 1 + h (l) + n ^ (3) 


Since the kernel functions are positive, the density calculated from the assumed interpolation 


has a lower bound The kernel functions assumed here are 
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where /i^ is the effective radius of the ith particle, is the particle separation distance, 


A denotes the unit step function, and a is a constant which allows for close packing of the 
particles at the reference bulk density. 

The hexahedral elements employed here are described elsewhere (e.g. Hallquist, 1983), 
so their interpolation is not discussed, except to note that hourglass deformation modes are 
avoided, via: (1) the use of subelement volumes to calculate hydrostatic tension, and (2) the 
use of a kernel based density function to describe compressed states. 


6 Kinetic and internal energy 


The Hamiltonian for the modeled system ( H ) is the sum of the kinetic ( T ) and internal ( U ) 
energy functions 


The kinetic energy is 


H = T + U 


T-± V>-ipW2 

2=1 1 


( 6 ) 

(7) 
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where the particle momenta are 

pt‘) = 77l W C (l) (8) 

with the velocity of the ith particle. The internal energy for the system is the sum 

U — Uq + U\ + U 2 ( 0 ) 


where the first term represents the stored energy in the particles. 

Up = y (10) 

t=i 

the second term represents stored energy due to hydrostatic tension in the subelements, 

U, = E £ U 1 - DM)Vi j ' k] K { Q j) (jW - l ) 2 A ( - 1) (11) 

j=u.=i 2 

and the third term represents stored energy due to elastic shear in the hexahedra. 

U 2 = E(1 - d 0) )Vo°Vo } tr[E e ^ T E e ^\ (12) 

with «0) the internal energy per unit mass for the particles, k,q ) and /x}/ ' bulk and shear 
moduli for the elements, and n e the number of elements. The assumed internal energy 
function has the general form 

U = U(c {i) .S il) .d u \D u) ,E p{]) ) (13) 


and defines the generalized conservative forces 



dU 
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(14) 


and the energy release rates 


G 0) = - 


dU 
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dU 
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dU 
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(15) 


for the Hamiltonian model, where 00) i s the thermodynamic temperature. Having defined 
the stored energy for the system, the next section describes the dissipative process models. 
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7 Plasticity, damage, and viscosity models 


Accurate hypervelocity impact models call for general descriptions of plastic flow, damage 
evolution, and material failure, and a numerical viscosity suitable for shock simulations. 
The plasticity model described here incorporates large strain kinematics, a rate dependent 
yield stress, and an isochoric plastic deformation constraint. The nonassociated flow rule is 


(Fahrenthold and Horban, 1997) 

EPb) = \U) {Itr[A wr A ( ^]}- 1/2 A w (16) 

with 

A ( d = C p(j) W (j) + W (j) C p(j) (17) 

where is a scalar multiplier and 

w( j ) = c p ^s^ + S U) C pU) _ -tr[C p(j) S {j) + S {j) C pU) ] I (18) 

3 

with the deviatoric stress tensor and plastic Cauchy-Green strain tensor defined by 

S ij) = (1 - d {j) ) 2 /i^E e0) , C pU) = 1 + 2 E p{j) (19) 

The assumed yield condition is 

jU) - T (i) _ Y {j) , t {]) = {itr[S 0)T S 0) ]} 1/2 (20) 

where is the effective stress. The yield stress is 

Y ij) = (i _ d (j) ){Yi j \ l + ^ ] e p(3) ) nU) +0[ ]) i U) }{\ - (21) 


where is the reference yield stress, e p ^ is the effective plastic strain, is the effective 
strain rate, 0^ is a strain hardening modulus, is a strain hardening exponent, 0^ is a 
strain rate hardening modulus, 0^ is a thermal softening modulus, and 6 H ^ is the maximum 
historical homologous temperature. The yield stress is limited to a specified maximum value 
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YML. The plastic strain increment at each time step is determined using a one step iteration 


procedure with 


AA (J) = 


( r (j) - Y {j) ) A [r {i) - Y {]) ) 


/ „\ \.~j 

(l-d0))2^ ) 

The second dissipative process modeled is evolution of the continuum damage. The 
simplest forms increment the damage over a fixed number of time steps (rif) once a failure 
criterion has been reached, for example 

A = — max{ K(max[a^\ t] — crp), A (max[e p ^\t\ — ip) } (23) 

n f 

where crb) is the maximum eigenvalue of the deviatoric stress tensor, a { p is a spall stress, 
eP is a failure strain, and the notation max\X. t ] denotes the maximum historical value of 
the argument X. Similarly for the normal damage 

A = — A (maxlP^K t] — P c ^) (24) 

n f 

where P ^ is the average tensile pressure in the element and P c(j ) is a critical pressure. 
Alternately normal damage evolution may be calculated using the rate dependent form 

Jj) 

D U) = -fr A (P 0) - P c(j) ) (25) 

h { p 

where cP is the local soundspeed and hp is the local particle radius, or 

A D {j) = (1 - D U) ) ^ Pi3) A (P (j) - P c{j) ) (26) 


where the damage increment at each time step serves to limit the tensile pressure, so as not 
to exceed the specified critical value. 

Finally a numerical viscosity is introduced. A variety of numerical viscosity models are 
discussed in the literature (Noh, 1978). Most take the general form 

fW = jp{Cp j) { c (i) - c 0) ) + Cp ]) \c (l) - c (j) |(c (t) - c (j) )} (27) 

where Cq'^ and C\ L ' 3} are numerical viscosity coefficients, nonzero only for neighboring par- 
ticles. 
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8 Entropy evolution 


The entropy evolution equations take the general form 

g(i) — Gjirr(i) _ gcon{i) (28) 

where the first term denotes irreversible entropy production due to the various dissipative 
processes and the second term is due to numerical heat conduction. These entropy evolution 
equations take the place of the energy balance relations included in conventional models of 
hypervelocity impact problems. The irreversible entropy production is calculated using 

n<*> -I 

^rr(t) = f (») . t (i) + ^2 - { + r i) (fc(0) + tr[G (k )T E pik ] ] } (29) 

where the are the element numbers associated with the ith particle, of number ny. The 
entropy flow due to numerical conduction is quantified using 

Q(i)gcon(i) _ - B {j) ) ( 30 ) 

j = i 

where i s a numerical conduction coefficient, nonzero only for neighboring particles. 

9 Hamilton’s equations 

The final form of Hamilton’s equations for the particle-element system may be derived using 
standard methods (Ginsberg, 1988), and consists of rate equations for the particle momentum 
and displacement states 

p(') = -g« - fW, c (l) = m^-y^ (31) 

augmented by the evolution equations for the particle entropies and the element damage and 
plastic variables. These are explicit equations, which may be integrated using Runge-Kutta 
or other methods to solve the impact problem. 
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10 ISC projectile simulations 


This section describes simulations of four different ISCL experiments, the latter performed 
by Grosch (1996 and 1997) on several different debris shield configurations. All of the 
simulations involved a projectile velocity of slightly over eleven kilometers per second, and 
all but one involved a velocity vector obliquity of 45 degrees. The ISCL projectiles were 
hollow aluminum cylinders with a length-to-diameter ratio less than two. and had a mass of 
approximately one gram. Since the projectile description was obtained from flash radiograph 
measurements, there is some uncertainty in the projectile mass and geometry data. 

In general the ISC projectiles exhibited both pitch and yaw with respect to the velocity 
vector, hence all of the simulations reported here are fully three dimensional. The models 
were composed of 100,000-500,000 particles and required as much as four days to simulate 
30-50 microseconds in physical time. The models were run in parallel on either 7 processors 
of an SGI Onyx or 32 processors of an SGI Origin, requiring up to 1GB of RAM. Computer 
resource constraints of course placed limits on the simulation times and the spatial resolutions 
of the models. Note that reducing the particle size by a factor of two would require a factor 
of eight increase in the number of particles and a factor of sixteen increase in the required 
wall clock time. Simulation parameters and material properties are listed in Appendices A 
and B. 

The first simulation involved a 45 degree oblique impact on an aluminum Whipple shield 
at a standoff distance of 7.62 cm. Figures la and lb show particle plots at impact and at 
46.6 microseconds after impact, while Figure lc shows an element plot of intact material 
at the simulation stop time. The simulation predicts a wall plate hole size (71 x 44 mm) 
somewhat greater than that observed in the experiment (60 x 20 mm). 

The second simulation involved a 45 degree oblique impact on an aluminum Whipple 
shield at a standoff distance of 11.43 cm. Figures 2a and 2b show particle plots at impact 
and at 45.0 microseconds after impact, while Figure 2c shows an element plot of intact 
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material at the simulation stop time. The simulation predicts a perforated region in the wall 
plate (25 x 10 mm) similar in size to the hole observed in the experiment (20 x 15 mm). 

The third simulation involved a normal impact on a dual plate aluminum shield at a 
standoff distance of 8.636 cm. In this case the axis of the cylindrical projectile and the 
velocity vector were significantly misaligned, again calling for a three dimensional simulation. 
Figures 3a and 3b show particle plots at impact and at 30.7 microseconds after impact, 
while Figure 3c shows an element plot of intact material at the simulation stop time. The 
simulation predicts a wall plate hole diameter (55 mm) somewhat greater than that observed 
in the experiment (44 mm). 

The fourth simulation involved a 45 degree oblique impact on a multilayer aluminum- 
Nextel-Kevlar shield at a standoff distance of 7.62 cm. Figures 4a and 4b show particle plots 
at impact and at 46.2 microseconds after impact, while Figure 4c shows an element plot 
of intact material at the simulation stop time. Consistent with the experimental results, 
the simulation predicts bulging but not perforation of the wall plate. It should be noted 
that some relevant material properties of Nextel and Kevlar are not well known, and are 
currently under study (Hiermaier et al., 1999). Although the linear elastic response of many 
composite materials has been well characterized, information on thermomechanical equation 
of state properties and plasticity properties is incomplete. Although the latter information 
is normally of secondary interest in structural design calculations, it is certainly of major 
interest in hypervelocity impact applications. 

The results just described indicate in general good agreement of the simulations with 
the experimental data. They do suggest a need for higher resolution models, longer physical 
simulation times, and better composite material models in future simulation work. 

11 Projectile shape effect 

As noted in the last section, the geometry of projectiles produced by ISCL experiments 
differs markedly from the solid spherical shape normally used in LGG tests. Since light 
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gas guns operate in a lower velocity regime, correlating the results of ISCL and LGG tests 
is complicated by an unknown projectile shape effect. In an attempt to investigate the 
significance of this projectile shape effect, the first three ISCL simulations described in the 
last section were repeated, with mass equivalent spherical projectiles replacing the actual 
hollow cylindrical ISC projectiles. Figures 5, 6, and 7 show element plots of the wall plate 
damage predictions obtained from simulations using hollow cylindrical ISCL projectiles and 
mass equivalent spherical projectiles, run in each case to the same simulation stop time. 

The results suggest that ISCL projectiles are more damaging than mass equivalent 
spheres, although the magnitude of the difference is difficult to quantify. In the first and 
third cases the projectile mass exceeds significantly the ballistic limit mass, and in all cases 
higher resolution models of the impact problems are needed in order to draw more definitive 
conclusions. However it should be noted that since the ISC projectiles: (1) are hollow, (2) 
exhibit pitch and yaw with respect to their velocity vector, and (3) involve rather low length- 
to-diameter ratios, one might expect to observe a modest projectile shape effect. Considering 
the complex nature of these highly oblique hypervelocity impact problems, it appears that 
more experimental and computational work is needed to address the question of projectile 
shape effects. 

12 Parallel speedup 

Three dimensional impact simulations require large memory and CPU time allocations. Pre- 
vious work on orbital debris shielding design (Faraud et al., 1999) has reported wall clock 
times as high as eighteen days for single processor simulations of three dimensional problems. 
Such turnaround times effectively preclude the use computer simulation in many engineering 
design projects. Parallel processing offers an opportunity to greatly reduce turnaround time 
and make three dimensional simulation a more practical design tool. 

The code used in the present work (Fahrenthold, 1999) was written for parallel execution 
on Onyx and Origin systems, using loop level compiler directives based on the OpenMP 
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standard. Alternative parallel implementations based on MPI constructs are in general more 
portable and presumably more efficient, although more difficult to implement. It should 
be emphasized that a high degree of parallelism must be present in the basic numerical 
algorithm, in order to achieve good speedup under any coding scheme. 

To evaluate parallel performance of the numerical algorithm and the code implementa- 
tion used here, speedup tests were run on Origin systems with up to 128 processors. The 
test problems were large (300,000 - 500,000 particles), to insure that a meaningful load was 
maintained on each CPU as the processor allocation increased. Figure 8 shows the absolute 
speedup measured for a 500,000 particle test problem, based on the wall clock time required 
for ten time steps at various CPU allocations. The dotted line shows the maximum theo- 
retical speedup, while the data points indicate the test results. At a CPU allocation of 64. 
the measured speedup is approximately two thirds of the theoretical maximum, indicating 
good parallel performance. At the maximum CPU allocation of 128. the efficiency drops to 
fifty percent. However the latter data point represents a factor of 64 reduction in wall clock 
time, indicating that a simulation which runs for over two months on one CPU can be run 
in one day on 128 processors. 

Massively parallel systems are characterized by distributed memory architectures, com- 
plicating somewhat the practical interpretation of speedup test data. The Origin system 
discussed here is composed of a collection of compute nodes , each of which consists of two 
processors and 512 MB of RAM. An individual user is allocated a discrete number of nodes 
for each particular job, that is allocations consisting of arbitrary combinations of processors 
and RAM are not permitted. As a result, a job which requires 1 GB of RAM will be allo- 
cated a minimum of four processors, and the meaningful speedup curve for such a problem 
is one measured relative to a CPU allocation of four. Figure 9 shows the results of a relative 
speedup test run on an Origin system, for ten time steps of a 300,000 particle test problem, 
using the code discussed in the present work. The solid line represents the maximum theo- 
retical relative speedup, while the data points show the test results. Again the data shows 
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good speedup for processor allocations as high as 64. 

High performance parallel computer systems are not yet commonplace in engineering 
design work. However the preceding results demonstrate that the numerical method used 
here can effectively exploit such resources, an important consideration as low cost, high 
performance parallel hardware becomes more widely available. 

13 Conclusion 

The present report describes a systematic test of the use of parallel computation and a 
hybrid particle-element algorithm to simulate a range of three dimensional orbital debris 
impact experiments. The numerical method appears to offer certain advantages in address- 
ing the three dimensional, multi-plate shield design problem. Additional work is needed 
to investigate model resolution, simulation time, projectile shape, and material property 
effects (including for example the use of the SESAME equation of state models). How- 
ever developments to date suggest that massively parallel computation using some type of 
mixed particle-continuum scheme offers excellent opportunities for significant advances in 
simulation-based debris shield design. 
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A Appendix: simulation parameters 

The simulated experiments are described in detail by Grosch (1996 and 1997). 


Simulation parameters 

Parameter 

Case 1 

Case 2 

Case 3 

Case 4 

SWRI Test Number 

7139-19 

7139-22 

7139-03 

7139-24 

Shield type 

Al Whipple 

Al Whipple 

Al dual plate 

Al-composite 

First aluminum plate thickness (cm) 

0.127 

0.127 

0.16002 

0.127 

Second aluminum plate thickness (cm) 

0.0 

0.0 

0.3175 

0.0 

Nextel areal density (g /cm 2 ) 

0.0 

0.0 

0.0 

0.4 

Kevlar areal density (g/ cm 2 ) 

0.0 

0.0 

0.0 

0.128 

Wall plate thickness (aluminum, cm) 

0.4826 

0.4826 

0.2032 

0.3175 

Maximum standoff (cm) 

7.62 

11.43 

8.636 

7.62 

Impact velocity (km/sec) 

11.41 

11.30 

11.16 

11.25 

Impact obliquity (velocity vector, deg) 

45 

45 

0 

45 

Projectile mass (aluminum, g) 

1.38 

0.85 

1.30 

1.07 

Projectile length-to-diameter ratio 

1.4 

1.2 

1.84 

1.1 

Projectile pitch (wrt velocity vector, deg) 

0 

11.6 

12.6 

0 

Projectile yaw (wrt velocity vector, deg) 

9.4 

19.3 

6.9 

0 

Number of particles 

142,867 

305,551 

265,251 

415,413 

Simulation time (fisec) 

46.6 

45.0 

30.7 

46.2 

Wall clock time (hours) 

28.1 

15.6 

15.6 

109.3 

Average number of processors 

6.9 

32 

32 

7.3 

System 

Onyx 

Origin 

Origin 

Onyx 


B Appendix: material properties 

Material properties were estimated using data from Steinberg (1996). Lee (1989), and Hiermaier et al. (1999). 
As indicated in the text, material models for the composites are the subject of current research. Listed below 
are the material properties used in the simulations. 


Material properties 

Parameter 

Aluminum 

Nextel 

Kevlar 

Equation of state type 

Mie-Gruneisen 

Linear 

Linear 

Shear modulus (Mbar) 

0.271 

0.164 

0.100 

Reference bulk density (g/cc) 

2.7 

0.82021 

0.741084 

Reference bulk modulus (Mbar) 

0.7832 

0.66633 

0.415389 

Initial yield stress (Mbar) 

0.0029 

0.008 

0.008 

Maximum yield stress (Mbar) 

0.0058 

0.008 

0.008 

Strain hardening exponent 

0.1 

0 

0 

Strain hardening modulus 

125.0 

0 

0 

Thermal softening modulus 

0.5 

1.0 

1.0 

Melt temperature (kilodegrees Kelvin) 

1.22 

1.22 1 

0.70 

Specific heat (Mbar-cm 3 per g-kilodegrees Kelvin) 

0.884 x 10~ 2 

0.884 x 10" 2 

1.420 x 10~ 2 

Spall stress (Mbar) 

0.012 

0.100 

0.100 

Plastic failure strain 

2.0 

0.2 

0.2 
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Figure la. Whipple shield with 7.62 cm standoff. 



Figure lc. Element plot at t = 46.6 microseconds. 
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Figure 2a. Whipple shield with 11.43 cm standoff. 



Figure 2b. Particle plot at t = 45.0 microseconds. 
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Figure 2c. Element plot at t = 45.0 microseconds. 



Figure 3a. Aluminum dual plate shield. 



Figure 3b. Particle plot at t = 30.7 microseconds. 
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Figure 3c. Element plot at t = 30.7 microseconds. 
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Figure 4a. Aluminum-Nextel-Kevlar shield. 



SWRI Test 7 1 39-24 (I - 46.2 microsec) 


Figure 4c. Element plot at t = 46.2 microseconds. 
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Figure 5a. Wall damage for ISC projectile, 
Whipple shield with 7.62 cm standoff. 
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Figure 5b. Wall damage for spherical projectile, 
Whipple shield with 7.62 cm standoff. 
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Figure 6a. Wall damage for ISC projectile, 
Whipple shield with 11.43 cm standoff. 
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Figure 6b. Wall damage for spherical projectile, 
Whipple shield with 11.43 cm standoff. 
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Figure 7a. Wall damage for ISC projectile, 
dual plate aluminum shield. 




Figure 7b. Wall damage for spherical projectile, 
dual plate aluminum shield. 
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Figure 9. Relative speedup for a problem with 300,000 particles. 
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